%% Description: CONCAVE dummy calculation
% Input: optiondata.mat file (cp_flag, strike price, implied volatility)
% Output: passdummy = CONCAVE indicator in the paper
%%
clear; clc;
%% load option data
load('optiondata.mat')
S=223.3;%stock price
rf=2e-4;%risk-free rate
ttm=8/365;%time-to-maturity
%% (1) Proceed only if there are at least 6 options with at least 2 calls and 2 puts
if height(optiondata)<6 || sum(strcmp(optiondata.cp_flag,'P'))<2 || sum(strcmp(optiondata.cp_flag,'C'))<2
    return
end
%% (2) Blend Implied Volatilites in the ATM region following Figlewski (2010)
atmboundL=S*0.98; atmboundH=S*1.02;  
locs=find(optiondata.strike_price>=atmboundL & optiondata.strike_price<=atmboundH);
strikesUlocs=unique(optiondata.strike_price(locs)); datanew=NaN(length(strikesUlocs),2); datanew(:,1)=strikesUlocs;
for strr=1:length(strikesUlocs)
    locsstrike=find(ismember(optiondata.strike_price,strikesUlocs(strr)));
    if length(locsstrike)==2
       wPut=(atmboundH-optiondata.strike_price(locsstrike(1)))./(atmboundH-atmboundL);
       putloc=find(ismember(optiondata.strike_price,strikesUlocs(strr))&ismember(optiondata.cp_flag,'P'));
       callloc=find(ismember(optiondata.strike_price,strikesUlocs(strr))&ismember(optiondata.cp_flag,'C'));
       datanew(strr,2)=wPut*optiondata.impl_volatility(putloc)+(1-wPut)*optiondata.impl_volatility(callloc);
       clear wPut
    end
end
optiondatav=optiondata;
optiondatav(locs,:)=[]; locnandatanew=find(isnan(datanew(:,2))); datanew(locnandatanew,:)=[];
fdata=[optiondatav.strike_price optiondatav.impl_volatility]; fdata=[datanew;fdata];
fdata=sortrows(fdata,1); fdata(:,1)=fdata(:,1)/S; % fdata has moneyness (K/S) on the first column and IVs on the second column
%% (3) Calculate the CONCAVE dummy variable (passdummy) 
% Loop through tolerance levels, start with a low RMSE and increase in small steps IF there is a negative IV point, a negative PDF point or more than 2 modes
 for fitloop=1:99
     TOLER=((.0001+(fitloop-1)*0.00005)^2)*size(fdata,1); clear peaks numpeaks prod2cn POINTS testarea domain
     %(a) interpolate IVs
     [iv,ki,order]=interp3c(fdata(:,1),fdata(:,2),TOLER);
     if sum(iv<0)>0
        continue
     end
     %(b) compute call option prices, the PDF and CDF following Breeden and Litzenberger (1978)
     calls=blsprice(S,ki.*S,rf,ttm,iv);
     [PDFemp,CDFemp,secder,firder]=deal(NaN(length(ki),1));
     for dens=2:length(ki)-1
         secder(dens)=(iv(dens+1)-2*iv(dens)+iv(dens-1))/((ki(dens+1)-ki(dens-1))^2);
         firder(dens)=(iv(dens+1)-iv(dens-1))/(ki(dens+1)-ki(dens-1));
         PDFemp(dens)=exp(rf*ttm)*(calls(dens+1)-2*calls(dens)+calls(dens-1))/((ki(dens+1)*S-ki(dens-1)*S)^2);
         CDFemp(dens)=exp(rf*ttm)*((calls(dens+1)-calls(dens-1))/(ki(dens+1)*S-ki(dens-1)*S))+1;
     end
     ki([1 1001])=[]; iv([1 1001])=[]; PDFemp([1 1001])=[]; CDFemp([1 1001])=[]; secder([1 1001])=[]; firder([1 1001])=[];
     peaks=findpeaks(PDFemp); numpeaks=length(peaks); locK2=find(abs(ki-fdata(2,1))==min(abs(ki-fdata(2,1)))); locK2=locK2(1);
     locKe1=find(abs(ki-fdata(end-1,1))==min(abs(ki-fdata(end-1,1)))); locKe1=locKe1(end);
     %(c) If there are only positive IVs, positive PDF points and number of modes<=2 break the tolerance loop and proceed           
     if sum(PDFemp<0)==0 && numpeaks<=2   
     %(d) Find concave regions, can be 0,1 or more and
     %set passdummy=1 (concave dummy=1) if the stationarity condition and minimum moneyness region of at least 0.03 are satisfied
     negsec=find(secder<0); diffsd=diff(negsec); locsbreaks=find(diffsd~=1); 
     if isempty(negsec)
        break
     elseif ~isempty(negsec) && isempty(locsbreaks)
        domain=ki(negsec(end))-ki(negsec(1));
        for fds=2:length(negsec)
            prod2cn(fds)=prod(firder(negsec(fds-1:fds)));
        end
        negprod=find(prod2cn<0); MINpoint=find(abs(firder(negsec))==min(abs(firder(negsec))));
        if domain>=0.03 && sum(prod2cn<0)>0 && ~isempty(intersect(negsec(negprod),locK2:locKe1)) || domain>=0.03 && ~isempty(intersect(negsec(MINpoint),locK2:locKe1)) && abs(firder(negsec(MINpoint)))<0.1
            passdummy=1;
            break
        else
            break
        end
      elseif ~isempty(negsec) && length(locsbreaks)==1
             domain(1)=ki(negsec(locsbreaks(1)))-ki(negsec(1)); domain(2)=ki(negsec(end))-ki(negsec(locsbreaks(1)+1)); 
             for fds=2:length(negsec)
                 prod2cn(fds)=prod(firder(negsec(fds-1:fds)));
             end
             prod2cn(locsbreaks(1)+1)=NaN; above15=find(domain>=0.03); negprod=find(prod2cn<0);
             if sum(domain>=0.03)>0
                if isequal(above15,1)
                   testarea=1:locsbreaks(1); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                   MINpoint=negsec(testarea(MINpoint));
                elseif isequal(above15,2)
                   testarea=locsbreaks(1)+1:length(negsec); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                   MINpoint=negsec(testarea(MINpoint));
                elseif isequal(above15,[1 2])
                   testarea=[1:locsbreaks(1) locsbreaks(1)+1:length(negsec)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:end);
                   MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                   MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                end            
                MINpoint=intersect(MINpoint,locK2:locKe1);
             end
             if sum(domain>=0.03)>0 && ~isempty(intersect(negsec(negprod),intersect(negsec(testarea),locK2:locKe1))) || sum(domain>=0.03)>0  && sum(abs(firder(MINpoint))<0.1)>0 
                passdummy=1;
                break
             else
                break
             end
     elseif ~isempty(negsec) && length(locsbreaks)==2
         domain(1)=ki(negsec(locsbreaks(1)))-ki(negsec(1)); domain(2)= ki(negsec(locsbreaks(2)))-ki(negsec(locsbreaks(1)+1));
         domain(3)=ki(negsec(end))-ki(negsec(locsbreaks(2)+1));
         for fds=2:length(negsec)
             prod2cn(fds)=prod(firder(negsec(fds-1:fds)));
         end
         prod2cn(locsbreaks(1)+1)=NaN; prod2cn(locsbreaks(2)+1)=NaN;  above15=find(domain>=0.03); negprod=find(prod2cn<0);
         if sum(domain>=0.03)>0
             if isequal(above15,1)
                 testarea=1:locsbreaks(1); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                 MINpoint=negsec(testarea(MINpoint));
             elseif isequal(above15,2)
                 testarea=locsbreaks(1)+1:locsbreaks(2); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                 MINpoint=negsec(testarea(MINpoint));
             elseif isequal(above15,3)
                 testarea=locsbreaks(2)+1:length(negsec);  MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                 MINpoint=negsec(testarea(MINpoint));
             elseif isequal(above15,[1 2])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[1 3])
                 testarea=[1:locsbreaks(1) locsbreaks(2)+1:length(negsec)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(2)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[2 3])
                 testarea=[locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:length(negsec)]; area1=negsec(locsbreaks(1)+1:locsbreaks(2)); area2=negsec(locsbreaks(2)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[1 2 3])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:length(negsec)]; area1=negsec(1:locsbreaks(1));
                 area2=negsec(locsbreaks(1)+1:locsbreaks(2)); area3=negsec(locsbreaks(2)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             end
             MINpoint=intersect(MINpoint,locK2:locKe1);
         end
         if  sum(domain>=0.03)>0 && ~isempty(intersect(negsec(negprod),intersect(negsec(testarea),locK2:locKe1))) || sum(domain>=0.03)>0  && sum(abs(firder(MINpoint))<0.1)>0
             passdummy=1;
             break
         else
             break
         end
     elseif ~isempty(negsec) && length(locsbreaks)==3
         domain(1)=ki(negsec(locsbreaks(1)))-ki(negsec(1)); domain(2)=ki(negsec(locsbreaks(2)))-ki(negsec(locsbreaks(1)+1));
         domain(3)=ki(negsec(locsbreaks(3)))-ki(negsec(locsbreaks(2)+1)); domain(4)=ki(negsec(end))-ki(negsec(locsbreaks(3)+1));
         for fds=2:length(negsec)
             prod2cn(fds)=prod(firder(negsec(fds-1:fds)));
         end
         prod2cn(locsbreaks(1)+1)=NaN; prod2cn(locsbreaks(2)+1)=NaN; prod2cn(locsbreaks(3)+1)=NaN; above15=find(domain>=0.03); negprod=find(prod2cn<0);
         if sum(domain>=0.03)>0
             if isequal(above15,1)
                 testarea=1:locsbreaks(1); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                 MINpoint=negsec(testarea(MINpoint));
             elseif isequal(above15,2)
                 testarea=locsbreaks(1)+1:locsbreaks(2); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                 MINpoint=negsec(testarea(MINpoint));
             elseif isequal(above15,3)
                 testarea=locsbreaks(2)+1:locsbreaks(3); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                 MINpoint=negsec(testarea(MINpoint));
             elseif isequal(above15,4)
                 testarea=locsbreaks(3)+1:length(negsec); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                 MINpoint=negsec(testarea(MINpoint));
             elseif isequal(above15,[1 2])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2)];  area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[1 3])
                 testarea=[1:locsbreaks(1) locsbreaks(2)+1:locsbreaks(3)];  area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(2)+1:locsbreaks(3));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[1 4])
                 testarea=[1:locsbreaks(1) locsbreaks(3)+1:length(negsec)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(3)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[2 3])
                 testarea=[locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:locsbreaks(3)]; area1=negsec(locsbreaks(1)+1:locsbreaks(2)); area2=negsec(locsbreaks(2)+1:locsbreaks(3));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[2 4])
                 testarea=[locsbreaks(1)+1:locsbreaks(2) locsbreaks(3)+1:length(negsec)];  area1=negsec(locsbreaks(1)+1:locsbreaks(2)); area2=negsec(locsbreaks(3)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[3 4])
                 testarea=[locsbreaks(2)+1:locsbreaks(3) locsbreaks(3)+1:length(negsec)]; area1=negsec(locsbreaks(2)+1:locsbreaks(3)); area2=negsec(locsbreaks(3)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[1 2 3])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:locsbreaks(3)];  area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2)); area3=negsec(locsbreaks(2)+1:locsbreaks(3));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[1 2 4])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2) locsbreaks(3)+1:length(negsec)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2)); area3=negsec(locsbreaks(3)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[1 3 4])
                 testarea=[1:locsbreaks(1) locsbreaks(2)+1:locsbreaks(3) locsbreaks(3)+1:length(negsec)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(2)+1:locsbreaks(3)); area3=negsec(locsbreaks(3)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[2 3 4])
                 testarea=[locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:locsbreaks(3) locsbreaks(3)+1:length(negsec)]; area1=negsec(locsbreaks(1)+1:locsbreaks(2)); area2=negsec(locsbreaks(2)+1:locsbreaks(3)); area3=negsec(locsbreaks(3)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[1 2 3 4])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:locsbreaks(3) locsbreaks(3)+1:length(negsec)];
                 area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2)); area3=negsec(locsbreaks(2)+1:locsbreaks(3)); area4=negsec(locsbreaks(3)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
                 MINpoint(4)=find(abs(firder(area4))==min(abs(firder(area4)))); MINpoint(4)=area4(MINpoint(4));
             end
             MINpoint=intersect(MINpoint,locK2:locKe1);
         end
         if sum(domain>=0.03)>0 && ~isempty(intersect(negsec(negprod),intersect(negsec(testarea),locK2:locKe1))) || sum(domain>=0.03)>0 && sum(abs(firder(MINpoint))<0.1)>0
             passdummy=1;
             break
         else
             break
         end
     elseif ~isempty(negsec) && length(locsbreaks)==4
         domain(1)=ki(negsec(locsbreaks(1)))-ki(negsec(1)); domain(2)=ki(negsec(locsbreaks(2)))-ki(negsec(locsbreaks(1)+1));
         domain(3)=ki(negsec(locsbreaks(3)))-ki(negsec(locsbreaks(2)+1)); domain(4)=ki(negsec(locsbreaks(4)))-ki(negsec(locsbreaks(3)+1));
         domain(5)=ki(negsec(end))-ki(negsec(locsbreaks(4)+1));
         for fds=2:length(negsec)
             prod2cn(fds)=prod(firder(negsec(fds-1:fds)));
         end
         prod2cn(locsbreaks(1)+1)=NaN; prod2cn(locsbreaks(2)+1)=NaN; prod2cn(locsbreaks(3)+1)=NaN; prod2cn(locsbreaks(4)+1)=NaN; above15=find(domain>=0.03); negprod=find(prod2cn<0);
         if sum(domain>=0.03)>0
             if isequal(above15,1)
                 testarea=1:locsbreaks(1); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                 MINpoint=negsec(testarea(MINpoint));
             elseif isequal(above15,2)
                 testarea=locsbreaks(1)+1:locsbreaks(2); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                 MINpoint=negsec(testarea(MINpoint));
             elseif isequal(above15,3)
                 testarea=locsbreaks(2)+1:locsbreaks(3); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                 MINpoint=negsec(testarea(MINpoint));
             elseif isequal(above15,4)
                 testarea=locsbreaks(3)+1:locsbreaks(4); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                 MINpoint=negsec(testarea(MINpoint));
             elseif isequal(above15,5)
                 testarea=locsbreaks(4)+1:length(negsec); MINpoint=find(abs(firder(negsec(testarea)))==min(abs(firder(negsec(testarea)))));
                 MINpoint=negsec(testarea(MINpoint));
             elseif isequal(above15,[1 2])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[1 3])
                 testarea=[1:locsbreaks(1) locsbreaks(2)+1:locsbreaks(3)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(2)+1:locsbreaks(3));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[1 4])
                 testarea=[1:locsbreaks(1) locsbreaks(3)+1:locsbreaks(4)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(3)+1:locsbreaks(4));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[1 5])
                 testarea=[1:locsbreaks(1) locsbreaks(4)+1:length(negsec)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[2 3])
                 testarea=[locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:locsbreaks(3)]; area1=negsec(locsbreaks(1):locsbreaks(2)); area2=negsec(locsbreaks(2)+1:locsbreaks(3));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[2 4])
                 testarea=[locsbreaks(1)+1:locsbreaks(2) locsbreaks(3)+1:locsbreaks(4)]; area1=negsec(locsbreaks(1):locsbreaks(2)); area2=negsec(locsbreaks(3)+1:locsbreaks(4));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[2 5])
                 testarea=[locsbreaks(1)+1:locsbreaks(2) locsbreaks(4)+1:length(negsec)];area1=negsec(locsbreaks(1):locsbreaks(2)); area2=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[3 4])
                 testarea=[locsbreaks(2)+1:locsbreaks(3) locsbreaks(3)+1:locsbreaks(4)]; area1=negsec(locsbreaks(2):locsbreaks(3)); area2=negsec(locsbreaks(3)+1:locsbreaks(4));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[3 5])
                 testarea=[locsbreaks(2)+1:locsbreaks(3) locsbreaks(4)+1:length(negsec)]; area1=negsec(locsbreaks(2):locsbreaks(3)); area2=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[4 5])
                 testarea=[locsbreaks(3)+1:locsbreaks(4) locsbreaks(4)+1:length(negsec)]; area1=negsec(locsbreaks(3):locsbreaks(4)); area2=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
             elseif isequal(above15,[1 2 3])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:locsbreaks(3)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2)); area3=negsec(locsbreaks(2)+1:locsbreaks(3));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[1 2 4])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2) locsbreaks(3)+1:locsbreaks(4)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2)); area3=negsec(locsbreaks(3)+1:locsbreaks(4));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[1 2 5])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2) locsbreaks(4)+1:length(negsec)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2)); area3=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[1 3 4])
                 testarea=[1:locsbreaks(1) locsbreaks(2)+1:locsbreaks(3) locsbreaks(3)+1:locsbreaks(4)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(2)+1:locsbreaks(3)); area3=negsec(locsbreaks(3)+1:locsbreaks(4));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[1 3 5])
                 testarea=[1:locsbreaks(1) locsbreaks(2)+1:locsbreaks(3) locsbreaks(4)+1:length(negsec)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(2)+1:locsbreaks(3)); area3=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[1 4 5])
                 testarea=[1:locsbreaks(1) locsbreaks(3)+1:locsbreaks(4) locsbreaks(4)+1:length(negsec)]; area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(3)+1:locsbreaks(4)); area3=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[2 3 4])
                 testarea=[locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:locsbreaks(3) locsbreaks(3)+1:locsbreaks(4)]; area1=negsec(locsbreaks(1)+1:locsbreaks(2)); area2=negsec(locsbreaks(2)+1:locsbreaks(3)); area3=negsec(locsbreaks(3)+1:locsbreaks(4));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[2 3 5])
                 testarea=[locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:locsbreaks(3) locsbreaks(4)+1:length(negsec)]; area1=negsec(locsbreaks(1)+1:locsbreaks(2)); area2=negsec(locsbreaks(2)+1:locsbreaks(3)); area3=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[2 4 5])
                 testarea=[locsbreaks(1)+1:locsbreaks(2) locsbreaks(3)+1:locsbreaks(4) locsbreaks(4)+1:length(negsec)]; area1=negsec(locsbreaks(1)+1:locsbreaks(2)); area2=negsec(locsbreaks(3)+1:locsbreaks(4)); area3=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[3 4 5])
                 testarea=[locsbreaks(2)+1:locsbreaks(3) locsbreaks(3)+1:locsbreaks(4) locsbreaks(4)+1:length(negsec)]; area1=negsec(locsbreaks(2)+1:locsbreaks(3)); area2=negsec(locsbreaks(3)+1:locsbreaks(4)); area3=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
             elseif isequal(above15,[1 2 3 4])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:locsbreaks(3) locsbreaks(3)+1:locsbreaks(4)];
                 area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2)); area3=negsec(locsbreaks(2)+1:locsbreaks(3)); area4=negsec(locsbreaks(3)+1:locsbreaks(4));
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
                 MINpoint(4)=find(abs(firder(area4))==min(abs(firder(area4)))); MINpoint(4)=area4(MINpoint(4));
             elseif isequal(above15,[1 2 3 5])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:locsbreaks(3) locsbreaks(4)+1:length(negsec)];
                 area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2)); area3=negsec(locsbreaks(2)+1:locsbreaks(3)); area4=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
                 MINpoint(4)=find(abs(firder(area4))==min(abs(firder(area4)))); MINpoint(4)=area4(MINpoint(4));
             elseif isequal(above15,[1 2 4 5])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2) locsbreaks(3)+1:locsbreaks(4) locsbreaks(4)+1:length(negsec)];
                 area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2)); area3=negsec(locsbreaks(3)+1:locsbreaks(4)); area4=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
                 MINpoint(4)=find(abs(firder(area4))==min(abs(firder(area4)))); MINpoint(4)=area4(MINpoint(4));
             elseif isequal(above15,[1 3 4 5])
                 testarea=[1:locsbreaks(1) locsbreaks(2)+1:locsbreaks(3) locsbreaks(3)+1:locsbreaks(4) locsbreaks(4)+1:length(negsec)];
                 area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(2)+1:locsbreaks(3)); area3=negsec(locsbreaks(3)+1:locsbreaks(4)); area4=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
                 MINpoint(4)=find(abs(firder(area4))==min(abs(firder(area4)))); MINpoint(4)=area4(MINpoint(4));
             elseif isequal(above15,[2 3 4 5])
                 testarea=[locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:locsbreaks(3) locsbreaks(3)+1:locsbreaks(4) locsbreaks(4)+1:length(negsec)];
                 area1=negsec(locsbreaks(1)+1:locsbreaks(2)); area2=negsec(locsbreaks(2)+1:locsbreaks(3)); area3=negsec(locsbreaks(3)+1:locsbreaks(4)); area4=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
                 MINpoint(4)=find(abs(firder(area4))==min(abs(firder(area4)))); MINpoint(4)=area4(MINpoint(4));
             elseif isequal(above15,[1 2 3 4 5])
                 testarea=[1:locsbreaks(1) locsbreaks(1)+1:locsbreaks(2) locsbreaks(2)+1:locsbreaks(3) locsbreaks(3)+1:locsbreaks(4) locsbreaks(4)+1:length(negsec)];
                 area1=negsec(1:locsbreaks(1)); area2=negsec(locsbreaks(1)+1:locsbreaks(2)); area3=negsec(locsbreaks(2)+1:locsbreaks(3)); area4=negsec(locsbreaks(3)+1:locsbreaks(4)); area5=negsec(locsbreaks(4)+1:end);
                 MINpoint(1)=find(abs(firder(area1))==min(abs(firder(area1)))); MINpoint(1)=area1(MINpoint(1));
                 MINpoint(2)=find(abs(firder(area2))==min(abs(firder(area2)))); MINpoint(2)=area2(MINpoint(2));
                 MINpoint(3)=find(abs(firder(area3))==min(abs(firder(area3)))); MINpoint(3)=area3(MINpoint(3));
                 MINpoint(4)=find(abs(firder(area4))==min(abs(firder(area4)))); MINpoint(4)=area4(MINpoint(4));
                 MINpoint(5)=find(abs(firder(area5))==min(abs(firder(area5)))); MINpoint(5)=area5(MINpoint(5));
             end
             MINpoint=intersect(MINpoint,locK2:locKe1);
         end
         if sum(domain>=0.03)>0 && ~isempty(intersect(negsec(negprod),intersect(negsec(testarea),locK2:locKe1))) || sum(domain>=0.03)>0  && sum(abs(firder(MINpoint))<0.1)>0
             passdummy=1;
             break
         else
             break
         end
     end
     clear iv ki otmcalls calls otmputs puts order negsec diffsd locsbreaks domain
     end
 end
